Data analysis for:
The Foraging-mode Paradigm: A Historical Overview Including a Reevaluation of its Predictions
Dylan Padilla, School of Life Sciences, Arizona State University, Tempe, AZ 85287, United States.
Library
library(AICcmodavg)
library(emo)
library(icons)
library(magick)
library(knitr)
library(MuMIn)
library(png)
library(raster)
library(scales)
library(vioplot)
library(xtable)
Loading data
data <- read.csv("data.csv")
head(data)> initial_mass final_mass strain env adult comments
> 1 0.0007245 NA s p NA died
> 2 0.0007720 NA s p NA died
> 3 0.0007095 NA s u NA died
> 4 0.0005566 NA s u NA died
> 5 0.0008130 NA s u NA died
> 6 0.0004055 0.0008495 s p NA <NA>data$growth <- data$final_mass - data$initial_mass
names(data)> [1] "initial_mass" "final_mass" "strain" "env" "adult" "comments" "growth"data <- data[ , c(3, 4, 7)]
data <- data[!is.na(data$growth), ]
data <- data[data$growth > 0, ]
table(data$strain, data$env)>
> p u
> r 29 29
> s 32 22data$strain[data$strain == "s"] <- "Sitter"
data$strain[data$strain == "r"] <- "Rover"
data$env[data$env == "p"] <- "Patchy"
data$env[data$env == "u"] <- "Lumpy"
data <- data[data$growth < 0.002, ] # Removing abnormal data
## Movement data
vid.dat <- read.csv('./OneDrive-2024-08-19/Results/Results_by_ind.csv', row.names = 1, sep = ';')
str(vid.dat)> 'data.frame': 92 obs. of 25 variables:
> $ Arena : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Sequence : chr "General" "General" "General" "General" ...
> $ Individual : chr "Ind0" "Ind0" "Ind0" "Ind0" ...
> $ Beginning_seq : num 0 0 0 0 0 0 0 0 0 0 ...
> $ End_seq : num 450 450 450 450 450 ...
> $ Prop_time_lost : num 0 0 0 0 0 0 0 0 0 0 ...
> $ Smoothing_filter_window : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Smoothing_Polyorder : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Moving_threshold : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Prop_time_moving : num 0.976 0.949 0.919 0.993 0.945 ...
> $ Average_Speed : num 0.518 2.36 5.962 5.412 2.979 ...
> $ Average_Speed_Moving : num 0.531 2.488 6.486 5.448 3.152 ...
> $ Traveled_Dist : num 233 1062 2683 2435 1341 ...
> $ Meander : num 1423 412 1387 1112 2120 ...
> $ Traveled_Dist_Moving : num 233 1062 2683 2435 1341 ...
> $ Meander_moving : num 1423 412 1387 1112 2120 ...
> $ Exploration_absolute_value : num 65.2 1148.3 150.6 1059.7 504.2 ...
> $ Exploration_relative_value : num 0.00671 0.12435 0.01568 0.11397 0.04999 ...
> $ Exploration_method : chr "Modern" "Modern" "Modern" "Modern" ...
> $ Exploration_area : int 10 10 10 10 10 10 10 10 10 10 ...
> $ Exploration_aspect_param : logi NA NA NA NA NA NA ...
> $ Mean_nb_neighbours : logi NA NA NA NA NA NA ...
> $ Prop_time_with_at_least_one_neighbour: logi NA NA NA NA NA NA ...
> $ Mean_shortest_interind_distance : logi NA NA NA NA NA NA ...
> $ Mean_sum_interind_distances : logi NA NA NA NA NA NA ...vid.dat <- vid.dat[c(2, 4, 1, 3, 5:92), ]
vid.dat[1:8, 1:5]> Arena Sequence Individual Beginning_seq End_seq
> Video SP - Apr 14 2024, 12 14 18.avi(3) 0 General Ind0 0 449.97
> Video RP - Apr 14 2024, 12 14 18.avi(4) 0 General Ind0 0 449.97
> Video SU - Apr 14 2024, 12 14 18.avi(1) 0 General Ind0 0 449.97
> Video RU - Apr 14 2024, 12 14 18.avi(2) 0 General Ind0 0 449.97
> Video SP - Apr 14 2024, 14 03 52.avi(1) 0 General Ind0 0 449.97
> Video RP - Apr 14 2024, 14 03 52.avi(2) 0 General Ind0 0 449.97
> Video SU - Apr 14 2024, 14 03 52.avi(3) 0 General Ind0 0 449.97
> Video RU - Apr 14 2024, 14 03 52.avi(4) 0 General Ind0 0 449.97mov.dat <- read.csv('movement.csv')
str(mov.dat)> 'data.frame': 92 obs. of 6 variables:
> $ mass : num 0.001991 0.001548 0.000821 0.001388 0.001131 ...
> $ strain : chr "r" "r" "s" "s" ...
> $ env : chr "p" "u" "p" "u" ...
> $ dis_trav: logi NA NA NA NA NA NA ...
> $ video : int 1 1 1 1 2 2 2 2 3 3 ...
> $ comments: logi NA NA NA NA NA NA ...head(mov.dat)> mass strain env dis_trav video comments
> 1 0.0019910 r p NA 1 NA
> 2 0.0015480 r u NA 1 NA
> 3 0.0008210 s p NA 1 NA
> 4 0.0013880 s u NA 1 NA
> 5 0.0011315 r p NA 2 NA
> 6 0.0014595 r u NA 2 NAs <- mov.dat[mov.dat$strain == "s", ]
r <- mov.dat[mov.dat$strain == "r", ]
dim(s)> [1] 46 6dim(r)> [1] 46 6ord.dat <- data.frame()
for(i in 1:46){
idx <- s[i, ]
idx2 <- r[i, ]
obj <- rbind(idx, idx2)
ord.dat <- rbind(ord.dat, obj)
}
ord.dat[1:8, ]> mass strain env dis_trav video comments
> 3 0.0008210 s p NA 1 NA
> 1 0.0019910 r p NA 1 NA
> 4 0.0013880 s u NA 1 NA
> 2 0.0015480 r u NA 1 NA
> 7 0.0014475 s p NA 2 NA
> 5 0.0011315 r p NA 2 NA
> 8 0.0012095 s u NA 2 NA
> 6 0.0014595 r u NA 2 NAvid.dat[1:8, 1:5]> Arena Sequence Individual Beginning_seq End_seq
> Video SP - Apr 14 2024, 12 14 18.avi(3) 0 General Ind0 0 449.97
> Video RP - Apr 14 2024, 12 14 18.avi(4) 0 General Ind0 0 449.97
> Video SU - Apr 14 2024, 12 14 18.avi(1) 0 General Ind0 0 449.97
> Video RU - Apr 14 2024, 12 14 18.avi(2) 0 General Ind0 0 449.97
> Video SP - Apr 14 2024, 14 03 52.avi(1) 0 General Ind0 0 449.97
> Video RP - Apr 14 2024, 14 03 52.avi(2) 0 General Ind0 0 449.97
> Video SU - Apr 14 2024, 14 03 52.avi(3) 0 General Ind0 0 449.97
> Video RU - Apr 14 2024, 14 03 52.avi(4) 0 General Ind0 0 449.97dim(ord.dat)> [1] 92 6vid.dat$mass <- ord.dat$mass
vid.dat$strain <- ord.dat$strain
vid.dat$env <- ord.dat$env
## Proportion of area visited
area.dat <- read.csv('Proportion-of-area-visited.csv', sep = ';')
area.dat <- area.dat[c(2, 4, 1, 3, 5:92), ]
vid.dat$prop_area_visited <- area.dat[ , 2]
names(vid.dat)> [1] "Arena" "Sequence" "Individual"
> [4] "Beginning_seq" "End_seq" "Prop_time_lost"
> [7] "Smoothing_filter_window" "Smoothing_Polyorder" "Moving_threshold"
> [10] "Prop_time_moving" "Average_Speed" "Average_Speed_Moving"
> [13] "Traveled_Dist" "Meander" "Traveled_Dist_Moving"
> [16] "Meander_moving" "Exploration_absolute_value" "Exploration_relative_value"
> [19] "Exploration_method" "Exploration_area" "Exploration_aspect_param"
> [22] "Mean_nb_neighbours" "Prop_time_with_at_least_one_neighbour" "Mean_shortest_interind_distance"
> [25] "Mean_sum_interind_distances" "mass" "strain"
> [28] "env" "prop_area_visited"str(vid.dat)> 'data.frame': 92 obs. of 29 variables:
> $ Arena : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Sequence : chr "General" "General" "General" "General" ...
> $ Individual : chr "Ind0" "Ind0" "Ind0" "Ind0" ...
> $ Beginning_seq : num 0 0 0 0 0 0 0 0 0 0 ...
> $ End_seq : num 450 450 450 450 450 ...
> $ Prop_time_lost : num 0 0 0 0 0 0 0 0 0 0 ...
> $ Smoothing_filter_window : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Smoothing_Polyorder : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Moving_threshold : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Prop_time_moving : num 0.949 0.993 0.976 0.919 0.945 ...
> $ Average_Speed : num 2.36 5.412 0.518 5.962 2.979 ...
> $ Average_Speed_Moving : num 2.488 5.448 0.531 6.486 3.152 ...
> $ Traveled_Dist : num 1062 2435 233 2683 1341 ...
> $ Meander : num 412 1112 1423 1387 2120 ...
> $ Traveled_Dist_Moving : num 1062 2435 233 2683 1341 ...
> $ Meander_moving : num 412 1112 1423 1387 2120 ...
> $ Exploration_absolute_value : num 1148.3 1059.7 65.2 150.6 504.2 ...
> $ Exploration_relative_value : num 0.12435 0.11397 0.00671 0.01568 0.04999 ...
> $ Exploration_method : chr "Modern" "Modern" "Modern" "Modern" ...
> $ Exploration_area : int 10 10 10 10 10 10 10 10 10 10 ...
> $ Exploration_aspect_param : logi NA NA NA NA NA NA ...
> $ Mean_nb_neighbours : logi NA NA NA NA NA NA ...
> $ Prop_time_with_at_least_one_neighbour: logi NA NA NA NA NA NA ...
> $ Mean_shortest_interind_distance : logi NA NA NA NA NA NA ...
> $ Mean_sum_interind_distances : logi NA NA NA NA NA NA ...
> $ mass : num 0.000821 0.001991 0.001388 0.001548 0.001448 ...
> $ strain : chr "s" "r" "s" "r" ...
> $ env : chr "p" "p" "u" "u" ...
> $ prop_area_visited : num 0.124 0.114 0.007 0.016 0.05 0.147 0.003 0.041 0.24 0.17 ...data.frame(vid.dat[1:12, 10:11])> Prop_time_moving Average_Speed
> Video SP - Apr 14 2024, 12 14 18.avi(3) 0.9486629 2.3603480
> Video RP - Apr 14 2024, 12 14 18.avi(4) 0.9934810 5.4120290
> Video SU - Apr 14 2024, 12 14 18.avi(1) 0.9757760 0.5177664
> Video RU - Apr 14 2024, 12 14 18.avi(2) 0.9191051 5.9617525
> Video SP - Apr 14 2024, 14 03 52.avi(1) 0.9452552 2.9791838
> Video RP - Apr 14 2024, 14 03 52.avi(2) 0.9984443 1.2822174
> Video SU - Apr 14 2024, 14 03 52.avi(3) 0.9691088 0.3256189
> Video RU - Apr 14 2024, 14 03 52.avi(4) 0.9758501 5.5721335
> Video SP - Apr 14 2024, 15 13 46.avi(1) 0.9748870 2.7576846
> Video RP - Apr 14 2024, 15 13 46.avi(2) 0.9659234 1.3956371
> Video SU - Apr 14 2024, 15 13 46.avi(3) 0.9698496 2.8616386
> Video RU - Apr 14 2024, 15 13 46.avi(4) 0.9392548 3.1724128
Modeling foraging behavior and
growth
##vid.dat <- vid.dat[vid.dat$Prop_time_moving > 0.7, ] ## extreme value
vid.dat$env[vid.dat$env == 'u'] <- 'Lumpy'
vid.dat$env[vid.dat$env == 'p'] <- 'Patchy'
vid.dat$strain[vid.dat$strain == 'r'] <- 'Rover'
vid.dat$strain[vid.dat$strain == 's'] <- 'Sitter'
str(vid.dat)> 'data.frame': 92 obs. of 29 variables:
> $ Arena : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Sequence : chr "General" "General" "General" "General" ...
> $ Individual : chr "Ind0" "Ind0" "Ind0" "Ind0" ...
> $ Beginning_seq : num 0 0 0 0 0 0 0 0 0 0 ...
> $ End_seq : num 450 450 450 450 450 ...
> $ Prop_time_lost : num 0 0 0 0 0 0 0 0 0 0 ...
> $ Smoothing_filter_window : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Smoothing_Polyorder : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Moving_threshold : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Prop_time_moving : num 0.949 0.993 0.976 0.919 0.945 ...
> $ Average_Speed : num 2.36 5.412 0.518 5.962 2.979 ...
> $ Average_Speed_Moving : num 2.488 5.448 0.531 6.486 3.152 ...
> $ Traveled_Dist : num 1062 2435 233 2683 1341 ...
> $ Meander : num 412 1112 1423 1387 2120 ...
> $ Traveled_Dist_Moving : num 1062 2435 233 2683 1341 ...
> $ Meander_moving : num 412 1112 1423 1387 2120 ...
> $ Exploration_absolute_value : num 1148.3 1059.7 65.2 150.6 504.2 ...
> $ Exploration_relative_value : num 0.12435 0.11397 0.00671 0.01568 0.04999 ...
> $ Exploration_method : chr "Modern" "Modern" "Modern" "Modern" ...
> $ Exploration_area : int 10 10 10 10 10 10 10 10 10 10 ...
> $ Exploration_aspect_param : logi NA NA NA NA NA NA ...
> $ Mean_nb_neighbours : logi NA NA NA NA NA NA ...
> $ Prop_time_with_at_least_one_neighbour: logi NA NA NA NA NA NA ...
> $ Mean_shortest_interind_distance : logi NA NA NA NA NA NA ...
> $ Mean_sum_interind_distances : logi NA NA NA NA NA NA ...
> $ mass : num 0.000821 0.001991 0.001388 0.001548 0.001448 ...
> $ strain : chr "Sitter" "Rover" "Sitter" "Rover" ...
> $ env : chr "Patchy" "Patchy" "Lumpy" "Lumpy" ...
> $ prop_area_visited : num 0.124 0.114 0.007 0.016 0.05 0.147 0.003 0.041 0.24 0.17 ...## Proportion of area covered
mod <- lm(prop_area_visited ~ mass*strain*env, data = vid.dat)
summary(mod)>
> Call:
> lm(formula = prop_area_visited ~ mass * strain * env, data = vid.dat)
>
> Residuals:
> Min 1Q Median 3Q Max
> -0.15443 -0.06215 -0.02083 0.06943 0.21500
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 0.11723 0.07979 1.469 0.1455
> mass 18.15868 54.18088 0.335 0.7383
> strainSitter -0.20760 0.12371 -1.678 0.0970 .
> envPatchy 0.28803 0.13347 2.158 0.0338 *
> mass:strainSitter 95.17618 88.27340 1.078 0.2840
> mass:envPatchy -152.05488 95.00854 -1.600 0.1133
> strainSitter:envPatchy -0.04573 0.18349 -0.249 0.8038
> mass:strainSitter:envPatchy 6.03000 134.26488 0.045 0.9643
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 0.09032 on 84 degrees of freedom
> Multiple R-squared: 0.3546, Adjusted R-squared: 0.3008
> F-statistic: 6.592 on 7 and 84 DF, p-value: 3.273e-06mod1 <- lm(prop_area_visited ~ mass*strain + env, data = vid.dat)
summary(mod1)>
> Call:
> lm(formula = prop_area_visited ~ mass * strain + env, data = vid.dat)
>
> Residuals:
> Min 1Q Median 3Q Max
> -0.15481 -0.06221 -0.02851 0.07036 0.21089
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 0.19850 0.06548 3.031 0.003209 **
> mass -34.18583 44.90920 -0.761 0.448584
> strainSitter -0.19425 0.08952 -2.170 0.032746 *
> envPatchy 0.06584 0.01925 3.420 0.000956 ***
> mass:strainSitter 71.56717 65.26893 1.096 0.275887
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 0.09156 on 87 degrees of freedom
> Multiple R-squared: 0.3131, Adjusted R-squared: 0.2815
> F-statistic: 9.915 on 4 and 87 DF, p-value: 1.173e-06mod2 <- lm(prop_area_visited ~ strain*env, data = vid.dat)
summary(mod2)>
> Call:
> lm(formula = prop_area_visited ~ strain * env, data = vid.dat)
>
> Residuals:
> Min 1Q Median 3Q Max
> -0.13956 -0.06178 -0.03193 0.07027 0.21578
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 0.14322 0.01904 7.521 4.36e-11 ***
> strainSitter -0.08291 0.02693 -3.079 0.00277 **
> envPatchy 0.08135 0.02693 3.021 0.00330 **
> strainSitter:envPatchy -0.03087 0.03808 -0.811 0.41981
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 0.09132 on 88 degrees of freedom
> Multiple R-squared: 0.3088, Adjusted R-squared: 0.2852
> F-statistic: 13.1 on 3 and 88 DF, p-value: 3.767e-07anova(mod2)> Analysis of Variance Table
>
> Response: prop_area_visited
> Df Sum Sq Mean Sq F value Pr(>F)
> strain 1 0.22246 0.222463 26.674 1.483e-06 ***
> env 1 0.09992 0.099924 11.981 0.0008313 ***
> strain:env 1 0.00548 0.005479 0.657 0.4198067
> Residuals 88 0.73391 0.008340
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod3 <- lm(prop_area_visited ~ strain + env, data = vid.dat)
summary(mod3)>
> Call:
> lm(formula = prop_area_visited ~ strain + env, data = vid.dat)
>
> Residuals:
> Min 1Q Median 3Q Max
> -0.14694 -0.06400 -0.02426 0.07215 0.20806
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 0.15093 0.01646 9.170 1.67e-14 ***
> strainSitter -0.09835 0.01901 -5.175 1.40e-06 ***
> envPatchy 0.06591 0.01901 3.468 0.00081 ***
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 0.09115 on 89 degrees of freedom
> Multiple R-squared: 0.3036, Adjusted R-squared: 0.288
> F-statistic: 19.4 on 2 and 89 DF, p-value: 1.015e-07anova(mod3)> Analysis of Variance Table
>
> Response: prop_area_visited
> Df Sum Sq Mean Sq F value Pr(>F)
> strain 1 0.22246 0.222463 26.778 1.399e-06 ***
> env 1 0.09992 0.099924 12.028 0.0008096 ***
> Residuals 89 0.73939 0.008308
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Model selection procedure based on AIC criterion
output.models <- model.sel(mod1, mod2, mod3)
sel.table <- round(as.data.frame(output.models)[7:11], 3)
names(sel.table)[1] <- "K"
sel.table$Model <- rownames(sel.table)
sel.table <- sel.table[ , c(6, 1, 2, 3, 4, 5)]
sel.table> Model K logLik AICc delta weight
> mod3 mod3 4 91.349 -174.237 0.000 0.604
> mod2 mod2 5 91.691 -172.684 1.554 0.278
> mod1 mod1 6 91.980 -170.972 3.265 0.118xtable(sel.table, digits = 3, caption = 'Caption here')> % latex table generated in R 4.4.1 by xtable 1.8-4 package
> % Sat Aug 24 10:49:20 2024
> \begin{table}[ht]
> \centering
> \begin{tabular}{rlrrrrr}
> \hline
> & Model & K & logLik & AICc & delta & weight \\
> \hline
> mod3 & mod3 & 4.000 & 91.349 & -174.237 & 0.000 & 0.604 \\
> mod2 & mod2 & 5.000 & 91.691 & -172.684 & 1.554 & 0.278 \\
> mod1 & mod1 & 6.000 & 91.980 & -170.972 & 3.265 & 0.118 \\
> \hline
> \end{tabular}
> \caption{Caption here}
> \end{table}## Model diagnosis
layout(matrix(c(0, 0, 0, 0,
1, 1, 2, 2,
1, 1, 2, 2,
0, 0, 0, 0), nrow = 4, ncol = 4, byrow = TRUE))
## Checking homogeneity of variance
plot(fitted(mod3), resid(mod3), col = "grey", pch = 20, xlab = "Fitted", ylab = "Residual", main = "Fitted versus Residuals", type = 'n', las = 1)
grid()
par(new = TRUE)
plot(fitted(mod3), resid(mod3), col = "grey", pch = 20, xlab = "Fitted", ylab = "Residual", main = "Fitted versus Residuals", las = 1)
abline(h = 0, col = "darkorange", lwd = 2)
## Checking normality
qqnorm(resid(mod3), col = "darkgrey", type = 'n', las = 1)
grid()
par(new = TRUE)
qqnorm(resid(mod3), col = "darkgrey", las = 1)
qqline(resid(mod3), col = "dodgerblue", lwd = 2)
## Distance traveled
mod4 <- lm(Traveled_Dist_Moving ~ mass*strain + env, data = vid.dat)
summary(mod4)>
> Call:
> lm(formula = Traveled_Dist_Moving ~ mass * strain + env, data = vid.dat)
>
> Residuals:
> Min 1Q Median 3Q Max
> -1009.14 -424.53 -53.04 350.55 1555.33
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 566.2 396.1 1.429 0.1565
> mass 362471.5 271648.8 1.334 0.1856
> strainSitter -169.6 541.5 -0.313 0.7549
> envPatchy 210.1 116.5 1.804 0.0747 .
> mass:strainSitter -287265.1 394801.6 -0.728 0.4688
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 553.8 on 87 degrees of freedom
> Multiple R-squared: 0.2534, Adjusted R-squared: 0.2191
> F-statistic: 7.381 on 4 and 87 DF, p-value: 3.629e-05mod5 <- lm(Traveled_Dist_Moving ~ strain*env, data = vid.dat)
summary(mod5)>
> Call:
> lm(formula = Traveled_Dist_Moving ~ strain * env, data = vid.dat)
>
> Residuals:
> Min 1Q Median 3Q Max
> -1130.4 -411.9 -105.2 299.3 1534.0
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 1148.55 115.07 9.981 3.93e-16 ***
> strainSitter -715.72 162.74 -4.398 3.05e-05 ***
> envPatchy 53.19 162.74 0.327 0.745
> strainSitter:envPatchy 278.99 230.15 1.212 0.229
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 551.9 on 88 degrees of freedom
> Multiple R-squared: 0.2501, Adjusted R-squared: 0.2245
> F-statistic: 9.783 on 3 and 88 DF, p-value: 1.233e-05anova(mod5)> Analysis of Variance Table
>
> Response: Traveled_Dist_Moving
> Df Sum Sq Mean Sq F value Pr(>F)
> strain 1 7636784 7636784 25.0742 2.813e-06 ***
> env 1 853966 853966 2.8039 0.09759 .
> strain:env 1 447556 447556 1.4695 0.22867
> Residuals 88 26801948 304568
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod6 <- lm(Traveled_Dist_Moving ~ strain + env, data = vid.dat)
summary(mod6)>
> Call:
> lm(formula = Traveled_Dist_Moving ~ strain + env, data = vid.dat)
>
> Residuals:
> Min 1Q Median 3Q Max
> -1060.7 -414.0 -109.2 366.2 1603.8
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 1078.80 99.92 10.797 < 2e-16 ***
> strainSitter -576.22 115.38 -4.994 2.92e-06 ***
> envPatchy 192.69 115.38 1.670 0.0984 .
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 553.3 on 89 degrees of freedom
> Multiple R-squared: 0.2376, Adjusted R-squared: 0.2204
> F-statistic: 13.87 on 2 and 89 DF, p-value: 5.727e-06anova(mod6)> Analysis of Variance Table
>
> Response: Traveled_Dist_Moving
> Df Sum Sq Mean Sq F value Pr(>F)
> strain 1 7636784 7636784 24.9426 2.921e-06 ***
> env 1 853966 853966 2.7891 0.09842 .
> Residuals 89 27249504 306174
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod7 <- lm(Traveled_Dist_Moving ~ strain, data = vid.dat)
summary(mod7)>
> Call:
> lm(formula = Traveled_Dist_Moving ~ strain, data = vid.dat)
>
> Residuals:
> Min 1Q Median 3Q Max
> -1157.02 -495.06 -59.71 351.98 1507.44
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 1175.15 82.39 14.263 < 2e-16 ***
> strainSitter -576.22 116.52 -4.945 3.51e-06 ***
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 558.8 on 90 degrees of freedom
> Multiple R-squared: 0.2137, Adjusted R-squared: 0.2049
> F-statistic: 24.46 on 1 and 90 DF, p-value: 3.507e-06anova(mod7)> Analysis of Variance Table
>
> Response: Traveled_Dist_Moving
> Df Sum Sq Mean Sq F value Pr(>F)
> strain 1 7636784 7636784 24.456 3.507e-06 ***
> Residuals 90 28103469 312261
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Model selection procedure based on AIC criterion
output.models <- model.sel(mod4, mod5, mod6, mod7)
sel.table <- round(as.data.frame(output.models)[7:11], 3)
names(sel.table)[1] <- "K"
sel.table$Model <- rownames(sel.table)
sel.table <- sel.table[ , c(6, 1, 2, 3, 4, 5)]
sel.table> Model K logLik AICc delta weight
> mod6 mod6 4 -710.085 1428.630 0.000 0.371
> mod7 mod7 3 -711.505 1429.282 0.652 0.268
> mod5 mod5 5 -709.323 1429.344 0.714 0.260
> mod4 mod4 6 -709.121 1431.231 2.600 0.101## Model diagnosis
layout(matrix(c(0, 0, 0, 0,
1, 1, 2, 2,
1, 1, 2, 2,
0, 0, 0, 0), nrow = 4, ncol = 4, byrow = TRUE))
## Checking homogeneity of variance
plot(fitted(mod6), resid(mod6), col = "grey", pch = 20, xlab = "Fitted", ylab = "Residual", main = "Fitted versus Residuals", type = 'n', las = 1)
grid()
par(new = TRUE)
plot(fitted(mod6), resid(mod6), col = "grey", pch = 20, xlab = "Fitted", ylab = "Residual", main = "Fitted versus Residuals", las = 1)
abline(h = 0, col = "darkorange", lwd = 2)
## Checking normality
qqnorm(resid(mod6), col = "darkgrey", type = 'n', las = 1)
grid()
par(new = TRUE)
qqnorm(resid(mod6), col = "darkgrey", las = 1)
qqline(resid(mod6), col = "dodgerblue", lwd = 2)
## Growth rate
int.model <- lm(growth ~ strain*env, data = data)
summary(int.model)>
> Call:
> lm(formula = growth ~ strain * env, data = data)
>
> Residuals:
> Min 1Q Median 3Q Max
> -3.238e-04 -1.061e-04 1.380e-06 1.033e-04 5.414e-04
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 3.221e-04 3.001e-05 10.735 < 2e-16 ***
> strainSitter 1.522e-05 4.569e-05 0.333 0.739690
> envPatchy -1.655e-04 4.321e-05 -3.830 0.000218 ***
> strainSitter:envPatchy 3.126e-05 6.242e-05 0.501 0.617629
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 0.0001616 on 105 degrees of freedom
> Multiple R-squared: 0.185, Adjusted R-squared: 0.1618
> F-statistic: 7.947 on 3 and 105 DF, p-value: 7.955e-05anova(int.model)> Analysis of Variance Table
>
> Response: growth
> Df Sum Sq Mean Sq F value Pr(>F)
> strain 1 7.4100e-09 7.4100e-09 0.2836 0.5955
> env 1 6.0859e-07 6.0859e-07 23.3073 4.712e-06 ***
> strain:env 1 6.5500e-09 6.5500e-09 0.2507 0.6176
> Residuals 105 2.7417e-06 2.6110e-08
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1noint.model <- lm(growth ~ strain + env, data = data)
summary(noint.model)>
> Call:
> lm(formula = growth ~ strain + env, data = data)
>
> Residuals:
> Min 1Q Median 3Q Max
> -0.0003334 -0.0001039 0.0000036 0.0001027 0.0005486
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 3.149e-04 2.622e-05 12.010 < 2e-16 ***
> strainSitter 3.196e-05 3.102e-05 1.030 0.305
> envPatchy -1.506e-04 3.107e-05 -4.845 4.35e-06 ***
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 0.000161 on 106 degrees of freedom
> Multiple R-squared: 0.1831, Adjusted R-squared: 0.1677
> F-statistic: 11.88 on 2 and 106 DF, p-value: 2.213e-05anova(noint.model)> Analysis of Variance Table
>
> Response: growth
> Df Sum Sq Mean Sq F value Pr(>F)
> strain 1 7.4100e-09 7.4100e-09 0.2856 0.5942
> env 1 6.0859e-07 6.0859e-07 23.4732 4.347e-06 ***
> Residuals 106 2.7482e-06 2.5930e-08
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Model selection procedure based on AIC criterion
output.models <- model.sel(int.model, noint.model)
sel.table <- round(as.data.frame(output.models)[5:9], 3)
names(sel.table)[1] <- "K"
sel.table$Model <- rownames(sel.table)
sel.table <- sel.table[ , c(6, 1, 2, 3, 4, 5)]
sel.table> Model K logLik AICc delta weight
> noint.model noint.model 4 798.862 -1589.340 0.000 0.725
> int.model int.model 5 798.992 -1587.402 1.938 0.275xtable(sel.table, digits = 3, caption = 'Caption here')> % latex table generated in R 4.4.1 by xtable 1.8-4 package
> % Sat Aug 24 10:49:20 2024
> \begin{table}[ht]
> \centering
> \begin{tabular}{rlrrrrr}
> \hline
> & Model & K & logLik & AICc & delta & weight \\
> \hline
> noint.model & noint.model & 4.000 & 798.862 & -1589.340 & 0.000 & 0.725 \\
> int.model & int.model & 5.000 & 798.992 & -1587.402 & 1.938 & 0.275 \\
> \hline
> \end{tabular}
> \caption{Caption here}
> \end{table}## Model diagnosis
layout(matrix(c(0, 0, 0, 0,
1, 1, 2, 2,
1, 1, 2, 2,
0, 0, 0, 0), nrow = 4, ncol = 4, byrow = TRUE))
## Checking homogeneity of variance
plot(fitted(noint.model), resid(noint.model), col = "grey", pch = 20, xlab = "Fitted", ylab = "Residual", main = "Fitted versus Residuals", type = 'n', las = 1)
grid()
par(new = TRUE)
plot(fitted(noint.model), resid(noint.model), col = "grey", pch = 20, xlab = "Fitted", ylab = "Residual", main = "Fitted versus Residuals", las = 1)
abline(h = 0, col = "darkorange", lwd = 2)
## Checking normality
qqnorm(resid(noint.model), col = "darkgrey", type = 'n', las = 1)
grid()
par(new = TRUE)
qqnorm(resid(noint.model), col = "darkgrey", las = 1)
qqline(resid(noint.model), col = "dodgerblue", lwd = 2)
## Proportion of area covered
xtable(data.frame(anova((mod3))), caption = 'Caption here...', digits = 3)> % latex table generated in R 4.4.1 by xtable 1.8-4 package
> % Sat Aug 24 10:49:21 2024
> \begin{table}[ht]
> \centering
> \begin{tabular}{rrrrrr}
> \hline
> & Df & Sum.Sq & Mean.Sq & F.value & Pr..F. \\
> \hline
> strain & 1 & 0.222 & 0.222 & 26.778 & 0.000 \\
> env & 1 & 0.100 & 0.100 & 12.028 & 0.001 \\
> Residuals & 89 & 0.739 & 0.008 & & \\
> \hline
> \end{tabular}
> \caption{Caption here...}
> \end{table}## Growth
xtable(formatC(as.matrix(anova(noint.model)), digits = 3, format = 'e'), caption = 'Caption here...', digits = 3)> % latex table generated in R 4.4.1 by xtable 1.8-4 package
> % Sat Aug 24 10:49:21 2024
> \begin{table}[ht]
> \centering
> \begin{tabular}{rlllll}
> \hline
> & Df & Sum Sq & Mean Sq & F value & Pr($>$F) \\
> \hline
> strain & 1.000e+00 & 7.405e-09 & 7.405e-09 & 2.856e-01 & 5.942e-01 \\
> env & 1.000e+00 & 6.086e-07 & 6.086e-07 & 2.347e+01 & 4.347e-06 \\
> Residuals & 1.060e+02 & 2.748e-06 & 2.593e-08 & NA & NA \\
> \hline
> \end{tabular}
> \caption{Caption here...}
> \end{table}
Plotting results of the model
## Barplot
## png('figure2.png', height = 7, width = 7, units = 'in', res = 350)
## Let's calculate the average value for each condition and each specie with the 'tapply' function
mean.area <- tapply(vid.dat$prop_area_visited, list(vid.dat$strain, vid.dat$env), mean)
as.matrix(mean.area)> Lumpy Patchy
> Rover 0.14321739 0.2245652
> Sitter 0.06030435 0.1107826## Plot boundaries
lim <- 2*max(mean.area)
## A function to add arrows on the chart
error.bar <- function(x, y, upper, lower = upper, length = 0.1,...){
arrows(x, y+upper, x, y-lower, angle = 90, code = 3, length = length, ...)
}
## Then I calculate the standard deviation for each specie and condition :
std.area <- tapply(vid.dat$prop_area_visited, list(vid.dat$strain, vid.dat$env), sd)
std.area> Lumpy Patchy
> Rover 0.10867171 0.07829944
> Sitter 0.07042819 0.10227012stdev.area <- as.matrix(std.area*1.96/46)
## I am ready to add the error bar on the plot using my "error bar" function!
cols <- c("black" , "gray")
par(mar = c(5, 6, 1, 1), mgp = c(3.8, 1, 0))
ze_barplot <- barplot(NA, beside = TRUE, legend.text = T, col = alpha(cols, 0.2), ylim = c(0, lim), ylab = expression(paste('Proportion of area visited')~(mm^2)), xlab = "Environment", las = 1)
grid()
par(new = TRUE)
ze_barplot <- barplot(mean.area, beside = TRUE, legend.text = T, col = alpha(cols, 0.5) , ylim = c(0, lim), ylab = expression(paste('Proportion of area visited')~(mm^2)), xlab = "Environment", las = 1)
error.bar(ze_barplot, mean.area, stdev.area, col = "black")
box()
## dev.off()## Barplot
## png('figure2.png', height = 7, width = 7, units = 'in', res = 350)
## Let's calculate the average value for each condition and each specie with the 'tapply' function
mean.gr <- tapply(data$growth, list(data$strain, data$env), mean)
as.matrix(mean.gr)> Lumpy Patchy
> Rover 0.0003221207 0.0001565889
> Sitter 0.0003373409 0.0002030645## Plot boundaries
lim <- 2*max(mean.gr)
## A function to add arrows on the chart
error.bar <- function(x, y, upper, lower = upper, length = 0.1,...){
arrows(x, y+upper, x, y-lower, angle = 90, code = 3, length = length, ...)
}
## Then I calculate the standard deviation for each specie and condition :
std.gr <- tapply(data$growth, list(data$strain, data$env), sd)
std.gr> Lumpy Patchy
> Rover 0.0001965322 9.225672e-05
> Sitter 0.0001966897 1.445087e-04stdev <- as.matrix(std.gr*1.96/46)
## I am ready to add the error bar on the plot using my "error bar" function!
cols <- c("black" , "gray")
par(mar = c(5, 6, 1, 1), mgp = c(3.8, 1, 0))
ze_barplot <- barplot(NA, beside = TRUE, legend.text = T, col = alpha(cols, 0.2), ylim = c(0, lim), ylab = expression(Growth~(Delta*gr%*%day^-1)), xlab = "Environment", las = 1)
grid()
par(new = TRUE)
ze_barplot <- barplot(mean.gr, beside = TRUE, legend.text = T, col = alpha(cols, 0.5) , ylim = c(0, lim), ylab = expression(Growth~(Delta*gr%*%day^-1)), xlab = "Environment", las = 1)
error.bar(ze_barplot, mean.gr, stdev, col = "black")
box()
## dev.off()
Figure 1
for(i in 1:23){
png(paste('figure', i, '.png', sep = ''), width = 7, height = 7, units = 'in', res = 350)
layout(matrix(c(1, 1, 2, 2,
1, 1, 2, 2,
3, 3, 4, 4,
3, 3, 4, 4), nrow = 4, ncol = 4, byrow = TRUE))
RP <- readPNG(source = paste('imgs/RP-path/figure', i, '.png', sep = ''))
SP <- readPNG(source = paste('imgs/SP-path/figure', i, '.png', sep = ''))
RU <- readPNG(source = paste('imgs/RU-path/figure', i, '.png', sep = ''))
SU <- readPNG(source = paste('imgs/SU-path/figure', i, '.png', sep = ''))
play <- readPNG('imgs/play-button.png')
par(mar = c(2.5, 1, 1, 0))
plot(NA, ylim = c(0, 1), xlim = c(0, 1), type = 'n', axes = FALSE, xlab = '', ylab = '', main = 'Rovers')
par(new = TRUE)
rasterImage(RP, 0, 0, 1, 1)
mtext('A', at = 0, line = -0.5)
par(mar = c(2.5, 1, 1, 0))
plot(NA, ylim = c(0, 1), xlim = c(0, 1), type = 'n', axes = FALSE, xlab = '', ylab = '', main = 'Sitters')
par(new = TRUE)
rasterImage(SP, 0, 0, 1, 1)
par(mar = c(2.5, 1, 1, 0))
plot(NA, ylim = c(0, 1), xlim = c(0, 1), type = 'n', axes = FALSE, xlab = '', ylab = '')
par(new = TRUE)
rasterImage(RU, 0, 0, 1, 1)
mtext('B', at = 0, line = -0.5)
par(mar = c(2.5, 1, 1, 0))
plot(NA, ylim = c(0, 1), xlim = c(0, 1), type = 'n', axes = FALSE, xlab = '', ylab = '')
par(new = TRUE)
rasterImage(SU, 0, 0, 1, 1)
par(xpd = TRUE)
rasterImage(play, -0.08, 0, 0.02, 0.11)
dev.off()
}
Figure 2
##png('figure2.png', width = 7, height = 7, units = 'in', res = 350)
layout(matrix(c(0, 1, 1, 0,
0, 1, 1, 0,
0, 2, 2, 0,
0, 2, 2, 0), nrow = 4, ncol = 4, byrow = TRUE))
par(mar = c(5, 6, 1, 1), mgp = c(3.8, 1, 0))
lim <- 2*max(mean.area)
ze_barplot <- barplot(NA, beside = TRUE, legend.text = T, col = alpha(cols, 0.2), ylim = c(0, lim), ylab = expression(paste('Proportion of area visited')~(mm^2)), xlab = "Environment", las = 1)
grid()
par(new = TRUE)
ze_barplot <- barplot(mean.area, beside = TRUE, legend.text = T, col = alpha(cols, 0.5) , ylim = c(0, lim), ylab = expression(paste('Proportion of area visited')~(mm^2)), xlab = "Environment", las = 1)
error.bar(ze_barplot, mean.area, stdev.area, col = "black")
box()
mtext('A', at = -0.45, line = -0.5)
par(mar = c(5, 6, 1, 1), mgp = c(3.8, 1, 0))
lim <- 2*max(mean.gr)
ze_barplot <- barplot(NA, beside = TRUE, legend.text = T, col = alpha(cols, 0.2), ylim = c(0, lim), ylab = expression(Growth~(Delta*gr%*%day^-1)), xlab = "Environment", las = 1)
grid()
par(new = TRUE)
ze_barplot <- barplot(mean.gr, beside = TRUE, legend.text = T, col = alpha(cols, 0.5) , ylim = c(0, lim), ylab = expression(Growth~(Delta*gr%*%day^-1)), xlab = "Environment", las = 1)
error.bar(ze_barplot, mean.gr, stdev, col = "black")
box()
mtext('B', at = -0.45, line = -0.5)
##dev.off()
Violin plots
##png("figure2-revised.png", height = 7, width = 7, units = "in", res = 360)
layout(matrix(c(0, 1, 1, 0,
0, 1, 1, 0,
0, 2, 2, 0,
0, 2, 2, 0), nrow = 4, ncol = 4, byrow = TRUE))
## Basic boxplot
par(mgp = c(2, 1, 0))
vioplot(prop_area_visited ~ strain + env, data = vid.dat, border = NA, method = "jitter", side = "right", ylab = expression(paste("Proportion of area visited")~(mm^2)), xlab = "Environment", col = "white", las = 1, xaxt = 'n', ylim = c(0, 0.45))
grid(nx = NULL, ny = NULL, col = alpha("lightgray", 0.5), lwd = 1, lty = 2)
par(new = TRUE)
vioplot(prop_area_visited ~ strain + env, data = vid.dat, border = NA, method = "jitter", side = "right", ylab = "", xlab = "", col = c(alpha("purple", 0.2), alpha("orange", 0.2)), las = 1, xaxt = 'n', ylim = c(0, 0.45))
axis(1, at = c(1.5, 3.5), labels = c('Uniform', 'Patchy'), cex = 0.8)
legend('topleft', legend = c('Rover', 'Sitter'), cex = 0.8, fill = c(alpha("purple", 0.2), alpha("orange", 0.2)), bty = 'n')
stripchart(prop_area_visited ~ strain + env, data = vid.dat, vertical = TRUE, method = "jitter", add = TRUE, pch = 20, col = c(alpha("purple", 0.3), alpha("orange", 0.5)), cex = 1.3)
mtext("A", side = 2, at = 0.5, line = 3, las = 1, font = 1)
## Basic boxplot
par(mgp = c(3, 1, 0))
vioplot(growth ~ strain + env, data = data, border = NA, method = "jitter", side = "right", ylab = expression(Growth~(Delta*gr%*%day^-1)), xlab = "Environment", col = "white", las = 1, xaxt = 'n')
grid(nx = NULL, ny = NULL, col = alpha("lightgray", 0.5), lwd = 1, lty = 2)
par(new = TRUE)
vioplot(growth ~ strain + env, data = data, border = NA, method = "jitter", side = "right", ylab = "", xlab = "", col = c(alpha("purple", 0.2), alpha("orange", 0.2)), las = 1, xaxt = 'n')
axis(1, at = c(1.5, 3.5), labels = c('Uniform', 'Patchy'), cex = 0.8)
stripchart(growth ~ strain + env, data = data, vertical = TRUE, method = "jitter", add = TRUE, pch = 20, col = c(alpha("purple", 0.3), alpha("orange", 0.5)), cex = 1.3)
mtext("B", side = 2, at = 0.00096, line = 3, las = 1, font = 1)
legend('topright', legend = c('Rover', 'Sitter'), cex = 0.8, fill = c(alpha("purple", 0.2), alpha("orange", 0.2)), bty = 'n')
##dev.off()